Домашнее задание №2 - Линейные модели. Градиентный спуск

В этом домашнем задании мы с вами научимся обучать линейные модели регрессии и классификации при помощи очень мощного, но в то же время довольно понятного алгоритма, который называется градиентный спуск. Помимо линейных моделей он используется и для обучения самых сложных нейронных сетей! Также мы потренируемся применять готовые реализации линейных моделей для задач регрессии и бинарной классификации.

import math
import sys
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import random
import pandas as pd

from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from sklearn.exceptions import NotFittedError
from sklearn.linear_model import (LinearRegression,
                                  LogisticRegression)
from sklearn.datasets import load_boston
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.metrics import (r2_score,
                             mean_squared_error)

from sklearn.datasets import make_classification
from sklearn.svm import SVC
from sklearn.metrics import precision_score, recall_score, f1_score, accuracy_score, RocCurveDisplay
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import plot_confusion_matrix

#!{sys.executable} -m pip install --upgrade scikit-learn

plt.rcParams["figure.figsize"] = 12, 9
sns.set_style("whitegrid")

SEED = 111
random.seed(SEED)
np.random.seed(SEED)

Маленькое теоретическое отступление

Основное свойство антиградиента (-1 * градиент) – он указывает в сторону наискорейшего убывания функции в данной точке. Соответственно, будет логично стартовать из некоторой точки, сдвинуться в сторону антиградиента, пересчитать антиградиент и снова сдвинуться в его сторону и т.д. Запишем это более формально.

Пусть $w_0$ – начальный набор параметров (коэффициентов линейной модели) ((например, нулевой или сгенерированный из некоторого, случайного распределения)). Тогда обычный градиентный спуск состоит в повторении следующих шагов до сходимости:

$$ w_{k + 1} = w_{k} - \eta \nabla_{w} Q(w_{k}), $$

где $\nabla_{w} Q(w_{k})$ – градиент функции потерь в точке $w_k$, а $\eta$ – скорость обучения (learning rate).

Градиентный спуск обычно останавливают, когда прошло заданное максимальное количество итераций или когда графиент близок к нулю (т.е. наши параметры практически не меняются). Для реализации второго варианта считают норму градиента (по сути длину вектора). Это можно сделать несколькими способами:

$$ l1_{norm} = \sum{|w_i|} $$$$ l2_{norm} = \sum{(w_i)^{2}} $$

Попробуем разобраться на простом примере. Рассмотрим функцию от двух переменных: $f(x, y) = \sin^2 x + \sin^2 y$

def f(w):
    """
    :param w: np.array(np.float) вектор из 2-х элементов
    :return: np.float
    """
    
    return np.sum(np.sin(w)**2)

Обратите внимание, что $x$ - numpy-array вектор длины 2.

Reminder:
Что мы хотим? Мы хотим найти минимум этой функции (в машинном обучении мы обычно хотим найти минимум функции потерь, например, MSE), а точнее найти $w_1$ и $w_2$ такие, что при них значение $f(w_1, w_2)$ минимально, то есть точку экстремума.

Как мы будем искать эту точку? Используем методы оптимизации (в нашем случае - минимизации). Одним из таких методов и является градиентный спуск.

Задание 1. Градиентный спуск для функции $f$ (1 балл)

Реализуйте функцию, которая будет осуществлять градиентный спуск для функции $f$:

Примечание: Вам нужно посчитать частные производные именно аналитически и переписать их в код, а не считать производные численно (через отношение приращения функции к приращению аргумента) -- в этих двух случаях могут различаться ответы, поэтому будьте внимательны.

def grad_f(w): 
    """
    Градиент функциии f, определенной выше.
        :param w: np.array[2]: float вектор из 2-х элементов
        :return: np.array[2]: float вектор из 2-х элементов
    """    
    # derivative from x
    derivative_x = 2 * np.sin(w[0])* np.cos(w[0])
    # derivative from y 
    derivative_y = 2 * np.sin(w[1])* np.cos(w[1])
    
    
    return np.array([derivative_x, derivative_y ])

Проверим, что градиент принимает вектор из двух чисел и выдает на этой точке верное значение

assert np.allclose(grad_f(np.array([1, 2])), 
                   np.array([0.90929743, -0.7568025])), "Что-то не так!"

def grad_descent_2d(f, grad_f, lr, num_iter=100, x0=None):
    """
    Функция, которая реализует градиентный спуск в минимум для функции f от двух переменных. 
        :param f: скалярная функция двух переменных
        :param grad_f: функция, возвращающая градиент функции f (устроена как реализованная вами выше grad_f)
        :param lr: learning rate алгоритма
        :param num_iter: количество итераций градиентного спуска
        :return: np.array[num_iter, 2] пар вида (x, f(x))
    """
    
    w0 = np.random.random(2)

    # будем сохранять значения аргументов и значений функции 
    # в процессе град. спуска в переменную history
    history = []

    # итерация цикла == шаг градиентнго спуска
    curr_w = w0.copy()
    for iter_num in range(num_iter):
        entry = np.hstack((curr_w, f(curr_w)))
        history.append(entry)
    
        curr_w -= lr * grad_f(curr_w)   #  YOUR CODE. Не забудьте про lr!

    return np.vstack(history)

Визуализируем точки градиентного спуска на 3D-графике нашей функции. Звездочками будут обозначены точки (тройки $w_1, w_2, f(w_1, w_2)$), по которым Ваш алгоритм градиентного спуска двигался к минимуму (Для того, чтобы написовать этот график, мы и сохраняли значения $cur\_w_1, cur\_w_2, f(cur\_w_1, cur\_w_2)$ в steps в процессе спуска).

Если у Вас правильно написана функция grad_descent_2d, то звездочки на картинке должны сходиться к одной из точек минимума функции. Вы можете менять начальные приближения алгоритма, значения lr и num_iter и получать разные результаты.

def gradient_desc_vis(f, grad_f, lr=0.01, num_iter=400):
    steps = grad_descent_2d(f, grad_f, lr=lr, num_iter=num_iter)
    
    X, Y = np.meshgrid(np.linspace(-3, 3, 100), np.linspace(-3, 3, 100))

    fig = plt.figure(figsize=(16, 10))
    ax = fig.gca(projection="3d")

    zs = np.array([f(np.array([x,y]))
                  for x, y in zip(np.ravel(X), np.ravel(Y))])
    Z = zs.reshape(X.shape)


    ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, zorder=2)

    ax.plot(xs=steps[:, 0], ys=steps[:, 1], zs=steps[:, 2],
            marker="*", markersize=5, zorder=3, 
            markerfacecolor="y", lw=3, c="black")

    ax.set_zlim(0, 5)
    ax.view_init(elev=60)
    plt.show()
    
    return steps

steps = gradient_desc_vis(f, grad_f)

<ipython-input-23-b28c0f7ff138>:22: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

We can see that the stars are indeed closer and closer to the mininmum!

Посмотрим на зависимость значения функции от шага градиентного спуска.

plt.figure(figsize=(14,7))
plt.xlabel("grad descent step number")
plt.ylabel("$f(x)$")
plt.title("Значение функции на каждом шаге гардиентного спуска.")

f_values = list(map(lambda x: x[2], steps))
plt.plot(f_values, label="gradient descent result")
plt.legend();

Задание 2. Реализация линейной регресии (суммарно 9 баллов)

Так как мы будем использовать градиентный спуск для обучения модели, важной частью является реализация функции потерь и функции для расчета ее градиента. Перем началом стоит напомнить, как считать градиент MSE. Вывод этой формулы можно найти здесь

$$ MSE = \frac{1}{N}\sum(y_{true} - y_{pred}) ^ 2 $$$$ \nabla{MSE} = \frac{2}{N} X^T (y_{pred} - y_{true}) $$

Здесь имеется в виду именно матричное умножение.

Задание 2.1. MSE и ее градиент (2 балла)

def mse(y_true, y_pred):
    """
    Функция потерь MSE.
        :param y_true: np.array[n_samples]: вектор из правильных ответов
        :param y_pred: np.array[n_samples]: вектор из предсказаний модели
        :return: значение функции потерь
    """
    
    if  y_true.shape[0] != y_pred.shape[0]:
        raise ValueError("Number of samples in both vectors should be equal")
        
    n = y_true.shape[0]
    vector_of_differences = y_true - y_pred
    sum_of_sq_differences = sum ([diff*diff for diff in vector_of_differences])
    mse = sum_of_sq_differences/n
    
    return mse


def mse_grad(y_true, y_pred, X):
    """
    Функция для расчета градиента MSE.
        :param y_true: np.array[n_samples]: вектор из правильных ответов
        :param y_pred: np.array[n_samples]: вектор из предсказаний модели
        :param X: np.array[n_samples, n_features]: матрица объекты x признаки
        :return: градиент функции потерь MSE
    """
    
    if  y_true.shape[0] != y_pred.shape[0]:
        raise ValueError("Number of samples in both vectors should be equal")
        
    n = y_true.shape[0]
    vector_of_differences = y_pred-y_true 
    x_transposed = X.transpose()
    grad = (2/n) * (x_transposed @ vector_of_differences)
    return grad


class MSELoss:
    """
    Класс, реализующий функцию потерь MSE. Нужен для того, чтобы
    объединять в одном месте функцию потерь и градиент для нее.
    """
    
    def __call__(self, y_true, y_pred):
        return mse(y_true, y_pred)
    
    def calculate_gradient(self, y_true, y_pred, X):
        return mse_grad(y_true, y_pred, X)

Мы будем использовать следующий класс для расчета градиента наших функций потерь:

class BasicGradientDescent:
    """
    Класс, позволяющий делать шаги градиентного спуска,
    а также рассчитывающих норму градиента.
    """
    
    def __init__(self, loss_function, grad_norm):
        self.loss = loss_function
        self.grad_norm = grad_norm
        
    
    def step(self, y, y_pred, X):
        grad_i = self.loss.calculate_gradient(y, y_pred, X)
        grad_i_norm = self._calculate_grad_norm(grad_i)
        
        return grad_i, grad_i_norm
            
            
    def _calculate_grad_norm(self, grad_i):
        if self.grad_norm == "l1":
            return np.abs(grad_i).sum()
        elif self.grad_norm == "l2":
            return np.sqrt(np.square(grad_i).sum())
        else:
            raise ValueError(f"I can't calculate {self.grad_norm} norm of gradient")

В данном задании нужно будет реализовать линейную регрессию и обучить ее при помощи градиентного спуска. Для этого нужно будет заполнять пропуски кода в соответствующих классах. Для начала мы реализуем базовый класс для всех линейных моделей, от которого потом будем наследоваться при реализации линейной и логистической регресий. Не переживайте, этот класс уже реализован, вам достостаточно просто разобраться с кодом.

class BaseLinearModel:
    """
    Класс, который представляет из себя базовую линейную модель, наследуюясь от которого, мы будем
    реализовывать линейную и логистическую регрессии.
    """
    
    def __init__(self, learning_rate, 
                 loss_function, fit_intercept,
                 n_iter, tol, optimizer, grad_norm):
        """
        Конструктор нашего класса. 
            :param learning_rate: скорость обучения
            :param loss_function: функция потерь (MSE или кросс-энтропия)
            :param fit_intercept: нужно ли нам включать свободных член в модель
            :param n_iter: количество итераций градиентного спуска
            :param tol: параметр для остановки градиентного спуска,
                        если норма градиента (l1 или l2) меньше tol, то останавливаемся
            :param optimizer: класс, который будет рассчитывать градиент и его норму
            :param grad_norm: тип нормы градиента l1 или l2
        """
        
        self.learning_rate = learning_rate
        self.loss = loss_function
        self.fit_intercept = fit_intercept
        self.n_iter = n_iter
        self.tol = tol
        self.grad_norm = grad_norm
        self.optimizer = optimizer(loss_function, grad_norm)
        
        # В начале параметры модели не заданы
        self.W = None
        
        
    def fit(self, X, y):
        """
        Метод для обучения нашей модели 
            :param X: матрица объекты x признаки
            :param y: вектор значений целевой переменной
            :return: обученная модель
        """
        
        # Сделаем из y вектор-столбец (n_samples, 1)
        y = y.reshape(-1, 1)
        n_samples = X.shape[0]
        
        # Добавим колонку из 1 в матрицу X
        if self.fit_intercept:
            ones_column = np.ones((n_samples, 1))
            X_new = np.hstack((ones_column, X))
        
        n_features = X_new.shape[1]
        
        # Инициализируем веса модели
        if self.W is None:
            self.W = np.random.randn(n_features, 1)
        
        # Обучаем модель градиентным спуском
        for i in range(self.n_iter):
            y_pred = self.predict(X)
            grad_i, grad_i_norm = self.optimizer.step(y, y_pred, X_new)
            
            # Если градиент близок к 0, останавливаемся
            if grad_i_norm <= self.tol:
                return self
            
            else:
                self.W -= self.learning_rate * grad_i
        
        return self
    
    
    def predict(self, X):
        raise NotImplementedError("It is a basic class for all linear models. You should implement it for descendant class.")
    
    
    def __repr__(self):
        return "Base linear model without prediction skill :("

Задание 2.2. Предсказания линейной регрессии (3 балла)

Реализуйте метод predict у класса CustomLinearRegression, не забудьте про свободный член!

class CustomLinearRegression(BaseLinearModel):
    def __init__(self, learning_rate: float = 1e-2, 
                 loss_function=MSELoss(), fit_intercept=True,
                 n_iter=1000, tol=1e-5, optimizer=BasicGradientDescent, grad_norm="l1"):
        
        # Если вы не проходили наследование и в частности `super`, то не страшно
        # коротко, с помощью этого мы можем вызывать методы родительского класса
        # в частности здесь мы используем метод `init`
        
        super().__init__(learning_rate=learning_rate, 
                         loss_function=loss_function, fit_intercept=fit_intercept,
                         n_iter=n_iter, tol=tol, optimizer=optimizer, grad_norm=grad_norm)
        
    
    def predict(self, X_test):
        """
        Метод для вычисления предсказаний 
            :param X_test: np.array[n_test_samples, n_features]: 
                           матрица объекты x признаки (тестовый датасет)
            :return: y_pred: np.array[n_test_samples, 1]: предсказания модели
        """
        
        if self.W is None:
            raise NotFittedError("This CustomLinearRegression instance is not fitted yet, run fit method.")
        
        n_test_samples = X_test.shape[0]
        
        if self.fit_intercept:
            ones_column = np.ones((n_test_samples, 1))
            X_test = np.hstack((ones_column, X_test))
            
        """
        YOUR CODE IS HERE
        """
        y_pred = X_test @ self.W
        return y_pred 
    
    
    def __repr__(self):
        return "My custom linear regression"

Проверим нашу реализацию на простом примере

X = np.random.randn(100, 1)
y = 2 * X + 5 + 0.5 * np.random.randn(100, 1)

plt.scatter(X, y);

custom_lin_reg = CustomLinearRegression()
custom_lin_reg.fit(X, y)

My custom linear regression

plt.scatter(X, y)
plt.plot(X, custom_lin_reg.predict(X));

Задание 2.3. Используем встроенную линейную регрессию (4 балла)

Поработаем с данными о ценах на дома в Бостоне. Постройте модель линейной регресии при помощи LinearRegression из sklearn. Не забудьте разделить данные на тренировочную и тестовую части, а также правильно предобработать признаки. В конце воспользуйтесь какими-то изученными метриками регресии и сделайте выводы о качестве полученной модели, а также о том, какие признаки наиболее важны с точки зрения полученной модели.

data = load_boston()
X, y = data["data"], data["target"]
feature_names = data["feature_names"]

# quick ugly look to check if data is standartized and all
f, axes = plt.subplots(figsize=(30,3), nrows=1, ncols=13)
for feature in range (X.shape[1]):
    axes_number = axes[feature]
    sns.histplot(x=X[:, feature], ax=axes_number)
    axes_number.set_title(data["feature_names"][feature])

Ваш ход:

x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
scaler = StandardScaler() # initialize a scaler 
# fit train data to scaler 
scaler.fit(x_train)
x_train_scaled = scaler.transform(x_train)
x_test_scaled = scaler.transform(x_test)

regr = LinearRegression()
regr.fit(x_train_scaled, y_train)
boston_y_pred = regr.predict(x_test_scaled)

for feature in range (X.shape[1]):
    print (data["feature_names"][feature], regr.coef_[feature])
    
from sklearn.metrics import (r2_score,
                             mean_squared_error)


plt.scatter(data["feature_names"], abs(regr.coef_), color="black")
plt.title('Regression coefficients absolute values for boston dataset', fontsize = '20')
plt.show()

CRIM -1.028001232211653
ZN 1.1638038505095065
INDUS -0.5261488458228758
CHAS 0.4429322755424035
NOX -1.472507368697788
RM 2.330718675029458
AGE 0.10091555043956957
DIS -3.047820682574309
RAD 2.7848307520849054
TAX -2.4820467440539282
PTRATIO -1.7900350206292446
B 0.8917244558137222
LSTAT -3.7513423750157084

<ipython-input-41-cc6fb4de38af>:14: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

The most important factors are DIS, LSTAT and RM. The least important - INDUS and AGE.

I tried to quickly visualize the dependencies of the predicted values from the x values (the picture looks tiny but can be easily enlarged by clicking twice). Here, we can see that RM or LSTAT indeed do look like they correlate with y variable!

f, axes = plt.subplots(figsize=(50,3), nrows=1, ncols=13)
for feature in range (X.shape[1]):
    axes_number = axes[feature]
    axes_number.scatter(x_test_scaled[:,feature], y_test, color="black", s=1)
    axes_number.scatter(x_test_scaled[:,feature], boston_y_pred, color="red", s=1)
    axes_number.set_title(data["feature_names"][feature])

# MSE - the smaller the better 
print(f"Mean squared error: {mean_squared_error(y_test, boston_y_pred)}")

# from - infinity to 1, 1 is the best 
print(f"Coefficient of determination: {r2_score(y_test, boston_y_pred)}")

Mean squared error: 30.1899905531444
Coefficient of determination: 0.6783233220896774

So, our R^2 is 0.68, which is not bad, but MSE is 30 - not exactly perfect.

Задание 3. Реализация логистической регресии (суммарно 10 баллов)

Логистическая регрессия не очень сильно отличается от обычной линейной регрессии и используется в задах классификации. Так как здесь мы снова будем пользоваться градиентным спуском, то нужно определить функцию потерь и ее градиент. Одним из самых популярных вариантов в задаче бинарной классификации является бинарная кросс-энтропия (BCE).

$$\mathcal L_{BCE}(y, \hat y) = -\sum_i \left[y_i\log\sigma(\hat y_i) + (1-y_i)\log(1-\sigma(\hat y_i))\right].$$

где $y$ это таргет желаемого результата и $\hat y$ является выходом модели. $\sigma$ - это логистическая функция, который преобразует действительное число $\mathbb R$ в вероятность $[0,1]$.

Единственная проблема данной функции это возможность получить 0 под знаком логарифма, что не очень хорошо. Попробуем справить с этим "в лоб". Скажем, что наши предсказания могут принимать значения от 0 + eps до 1 - eps, где eps очень маленькое число.

Задание 3.1. Реализация сигмоиды (0.5 баллов)

Реализуйте функцию sigmoid, которая переводит действительное число $\mathbb R$ в вероятность $[0,1]$.

def sigmoid(output):
    # output результат X@w (-inf, +inf)
    probability = 1 / (math.e ** (-output)+1)
    return probability

Задание 3.2. BCE Loss и ее градиент (2.5 балла)

Так как мы с вами только начинаем изучать машинное обучение, то было бы слишком жестоко просить вас вычислить градиент BCE Loss (он не так сложен, просто нужно привыкнуть). Поэтому сразу напишем формулу для него:

$$ \nabla{\mathcal L_{BCE}(y, \hat y), X} = X^T (\sigma({\hat{y}}) - y) $$

def bce(y_true, y_pred, eps=1e-15):
    """
    Функция потерь BCE.
        :param y_true: np.array[n_samples]: вектор из правильных ответов 0/1
        :param y_pred: np.array[n_samples]: вектор из предсказаний модели (вероятности)
        :return: значение функции потерь
    """
    
    if  y_true.shape[0] != y_pred.shape[0]:
        raise ValueError("Number of samples in both vectors should be equal")
        
    n = y_true.shape[0]
    
    # So I want escape log(0)
    y_pred = np.clip(y_pred, eps, 1 - eps)
    sum_bce = 0
    for i in range (n): 
        yi = y_true[i] 
        ypi = y_pred[i]
        sum_bce -= yi * math.log(ypi) + (1 - yi) * math.log(1 - ypi)
    bce = sum_bce/n
    return bce
    

def bce_grad(y_true, y_pred, X):
    """
    Функция потерь BCE.
        :param y_true: np.array[n_samples]: вектор из правильных ответов 0/1
        :param y_pred: np.array[n_samples]: вектор из предсказаний модели (вероятности)
        :param X: np.array[n_samples, n_features]: матрица объекты x признаки
        :return: значение функции потерь
    """
    
    if  y_true.shape[0] != y_pred.shape[0]:
        raise ValueError("Number of samples in both vectors should be equal")
        
    x_transposed = X.transpose()
    bce_grad = x_transposed @ (y_pred - y_true)
    return bce_grad


class BCELoss:
    """
    Класс, реализующий функцию потерь BCE. Нужен для того, чтобы
    объединять в одном месте функцию потерь и градиент для нее.
    """
    
    def __call__(self, y_true, y_pred):
        return bce(y_true, y_pred)
    
    def calculate_gradient(self, y_true, y_pred, X):
        return bce_grad(y_true, y_pred, X)

Задание 3.3. Предсказания логистической регрессии (2 балла)

Реализуйте метод predict у класса CustomLogisticRegression, не забудьте про свободный член!

class CustomLogisticRegression(BaseLinearModel):
    def __init__(self, learning_rate: float = 1e-3,
                 loss_function=BCELoss(), fit_intercept=True,
                 n_iter=1000, tol=1e-5, optimizer=BasicGradientDescent, grad_norm="l1"):
        
        super().__init__(learning_rate=learning_rate,
                         loss_function=loss_function, fit_intercept=fit_intercept,
                         n_iter=n_iter, tol=tol, optimizer=optimizer, grad_norm=grad_norm)
    
    
    def predict(self, X_test):
        if self.W is None:
            raise NotFittedError("This CustomLogisticRegression instance is not fitted, run fit method.")
        
        n_test_samples = X_test.shape[0]
        if self.fit_intercept:
            ones_column = np.ones((n_test_samples, 1))
            X_test = np.hstack((ones_column, X_test))
            
        """
        YOUR CODE IS HERE
        """
        y_pred = sigmoid(X_test @ self.W)
        return y_pred 
    
    
    def __repr__(self):
        return "My custom logistic regression"

Снова проверим работу алгоритма на простом примере

# Создадим датасет из 1 признака и 2 классов
X, y = make_classification(n_features=1, n_informative=1,
                           n_redundant=0, n_clusters_per_class=1)

plt.scatter(X, y);

custom_log_reg = CustomLogisticRegression()
custom_log_reg.fit(X, y)
y_pred = custom_log_reg.predict(X)

plt.scatter(X, y)
plt.scatter(X, y_pred);

Проверьте качество работы модели при помощи известных вам метрик бинарной классификации.

# turn the probability into a class, without the 'unknown' type
y_pred_classes = [1 if i > 0.5 else 0 for i in y_pred ]

# calculate accuracy and some other metrics 
print ('Accuracy: ', accuracy_score(y, y_pred_classes))
print ('Precision: ', precision_score(y,y_pred_classes))
print ('Recall: ', recall_score(y,y_pred_classes))
print ('F1: ', f1_score(y, y_pred_classes))

Accuracy:  0.98
Precision:  0.98
Recall:  0.98
F1:  0.98

Accuracy is 87 % - looks not bad!

Задание 3.4. Применение логистической регрессии (5 баллов)

Мы будем использовать данные по свойствам покемонов (https://www.kaggle.com/abcsds/pokemon). В данном задании вам необходимо сначала сделать краткий EDA (Посмотреть на данные и их распределения, а также посмотреть, как различные признаки связаны между собой и с целевой переменной (Legendary)).

pokemon = pd.read_csv("../data/Pokemon.csv")
pokemon.head()

	#	Name	Type 1	Type 2	Total	HP	Attack	Defense	Sp. Atk	Sp. Def	Speed	Generation	Legendary
0	1	Bulbasaur	Grass	Poison	318	45	49	49	65	65	45	1	False
1	2	Ivysaur	Grass	Poison	405	60	62	63	80	80	60	1	False
2	3	Venusaur	Grass	Poison	525	80	82	83	100	100	80	1	False
3	3	VenusaurMega Venusaur	Grass	Poison	625	80	100	123	122	120	80	1	False
4	4	Charmander	Fire	NaN	309	39	52	43	60	50	65	1	False

pokemon.describe().round(5)

	#	Total	HP	Attack	Defense	Sp. Atk	Sp. Def	Speed	Generation
count	800.00000	800.00000	800.00000	800.00000	800.0000	800.00000	800.00000	800.00000	800.00000
mean	362.81375	435.10250	69.25875	79.00125	73.8425	72.82000	71.90250	68.27750	3.32375
std	208.34380	119.96304	25.53467	32.45737	31.1835	32.72229	27.82892	29.06047	1.66129
min	1.00000	180.00000	1.00000	5.00000	5.0000	10.00000	20.00000	5.00000	1.00000
25%	184.75000	330.00000	50.00000	55.00000	50.0000	49.75000	50.00000	45.00000	2.00000
50%	364.50000	450.00000	65.00000	75.00000	70.0000	65.00000	70.00000	65.00000	3.00000
75%	539.25000	515.00000	80.00000	100.00000	90.0000	95.00000	90.00000	90.00000	5.00000
max	721.00000	780.00000	255.00000	190.00000	230.0000	194.00000	230.00000	180.00000	6.00000

We can see that data are not standartised, so we will need to do the standartisation.

Now we can look at the distribution of every feature, except for the name and number (they don't look exactly meaningful).

f, axes = plt.subplots(figsize=(60,8), nrows=1, ncols=11)
for feature in range (11):
    axes_number = axes[feature]
    sns.histplot(x=pokemon.iloc[:, feature+2].dropna(), ax=axes_number)

<__array_function__ internals>:5: RuntimeWarning: Converting input from bool to <class 'numpy.uint8'> for compatibility.
<__array_function__ internals>:5: RuntimeWarning: Converting input from bool to <class 'numpy.uint8'> for compatibility.

We can see that we have much less legendary guys than normal ones - should use it when splitting on train and test! But are features interconnectted somehow? We can visualize correlation matrix.

# correlation matrix
corr = pokemon.corr()
corr.style.background_gradient(cmap='coolwarm')

	#	Total	HP	Attack	Defense	Sp. Atk	Sp. Def	Speed	Generation	Legendary
#	1.000000	0.119813	0.097614	0.102298	0.094786	0.088759	0.085817	0.010733	0.982516	0.153396
Total	0.119813	1.000000	0.618748	0.736211	0.612787	0.747250	0.717609	0.575943	0.048384	0.501758
HP	0.097614	0.618748	1.000000	0.422386	0.239622	0.362380	0.378718	0.175952	0.058683	0.273620
Attack	0.102298	0.736211	0.422386	1.000000	0.438687	0.396362	0.263990	0.381240	0.051451	0.345408
Defense	0.094786	0.612787	0.239622	0.438687	1.000000	0.223549	0.510747	0.015227	0.042419	0.246377
Sp. Atk	0.088759	0.747250	0.362380	0.396362	0.223549	1.000000	0.506121	0.473018	0.036437	0.448907
Sp. Def	0.085817	0.717609	0.378718	0.263990	0.510747	0.506121	1.000000	0.259133	0.028486	0.363937
Speed	0.010733	0.575943	0.175952	0.381240	0.015227	0.473018	0.259133	1.000000	-0.023121	0.326715
Generation	0.982516	0.048384	0.058683	0.051451	0.042419	0.036437	0.028486	-0.023121	1.000000	0.079794
Legendary	0.153396	0.501758	0.273620	0.345408	0.246377	0.448907	0.363937	0.326715	0.079794	1.000000

It seems that Generation and index are interconnected, also total correlates with HP, speed of atack and speed of defence. We can look at the plots of how these variables loook like against each other in the pandas_profiling repoprt, section "Interactions"

import pandas_profiling as pp
pp.ProfileReport(pokemon) 

So let's use the dataset without the index number and name. We have to recode Type 1, Type 2 and legendary.

Pokemons can be of 18 different types. We can split both type 1 and type 2 into 18 different columns, but that will give us + 36 new features. It's likely they will not be very meaningful. Alternatively, we can turn each class into a number (from 0 to 18) - but that will infer that some classes are more alike than the others. In addition, almost half of the pokemons do not have Type 2 at all. I propose to get rid of the Type variable, recode Legendary in 1 - legendary and 0 - not legendary and standartize the train and test sets.

Мы будем предсказывать является ли покемон легендарным или нет. Замените логическое значение колонки на числовое (перекодировав на 0 и 1). Также подумайте, как в этом случае лучше закодировать категориальные признаки (может быть, лучше их просто выбросить?).

le = LabelEncoder()
pokemon.Legendary = le.fit_transform(pokemon.Legendary)
pokemon.head()
p_x = pokemon[['Total','HP', 'Attack', 'Defense', 'Sp. Atk', 'Sp. Def', 'Speed', 'Generation']]
p_y = pokemon[['Legendary']]
p_x.head()

	Total	HP	Attack	Defense	Sp. Atk	Sp. Def	Speed	Generation
0	318	45	49	49	65	65	45	1
1	405	60	62	63	80	80	60	1
2	525	80	82	83	100	100	80	1
3	625	80	100	123	122	120	80	1
4	309	39	52	43	60	50	65	1

Разделите ваши данные на тестовую и тренировочную выборку.

x_train, x_test, y_train, y_test = train_test_split(p_x, p_y , test_size=0.3, stratify=p_y)
scaler = StandardScaler() # initialize a scaler 
# fit train data to scaler 
scaler.fit(x_train)
x_train_scaled = scaler.transform(x_train)
x_test_scaled = scaler.transform(x_test)

Обучите модель LogisticRegression из sklearn.

regr = LogisticRegression()
regr.fit(x_train_scaled, y_train.values.ravel())
pokemon_y_pred = regr.predict(x_test_scaled)
pokemon_y_pred

array([0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
       0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0],
      dtype=int64)

Выведите метрики вашего классификатора:

1. Нарисуйте confusion matrix.

2. Изобразите ROC кривую и посчитайте площадь под ней.

3. Скажите, какие признаки оказались наиболее важны для модели.

Confusion matrix:

plot_confusion_matrix(regr, x_test_scaled, y_test)  
plt.show()

C:\Users\emeli\anaconda3\lib\site-packages\sklearn\utils\deprecation.py:87: FutureWarning: Function plot_confusion_matrix is deprecated; Function `plot_confusion_matrix` is deprecated in 1.0 and will be removed in 1.2. Use one of the class methods: ConfusionMatrixDisplay.from_predictions or ConfusionMatrixDisplay.from_estimator.
  warnings.warn(msg, category=FutureWarning)
<ipython-input-62-a00d0e37b14f>:2: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

Looks like most of the pokemons (216) are not legendary and predicted to be non-legendary!

ROC-curve

RocCurveDisplay.from_predictions(y_test, pokemon_y_pred)
plt.show()

<ipython-input-63-ce2705b589a8>:2: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

The square under the ROC-curve is 0.73 - ideally, would be one.

plt.scatter(p_x.columns, abs(regr.coef_), color="black")
plt.title('Regression coefficients absolute values for pokemon dataset', fontsize = '20')
plt.show()

<ipython-input-65-6953a6c75494>:3: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

Looks like Speed of attack, Speed and Total are the most important features!

Задание 4. Расскажите о вашей любимой музыкальной группе (исполнителе) (0.5 балла)

Расскажите, как вы познакомились с этой группой и скиньте несколько наиболее любимых треков)

Одна из моих любимых групп - "Немного нервно". Они исполняют дрим-фолк и прочую странную фиготу, а познакомила меня с этой группой лучшая подруга, когда мы обе были ещё в школе, а главная исполнительница ведёт твиттер под ником "Грушновато" (https://twitter.com/nervno_nemnogo). Примеры их песен - в видео ниже.

from IPython.display import YouTubeVideo

YouTubeVideo('sL-entmdnwA', width=800, height=500)

YouTubeVideo('_0Z3VmHmWKU', width=800, height=500)

Therapy time

Напишите здесь ваши впечатления о задании: было ли интересно, было ли слишком легко или наоборот сложно и тд. Также сюда можно написать свои идеи по улучшению заданий, а также предложить данные, на основе которых вы бы хотели построить следующие дз.

Ваши мысли:

Кажется, что второе задание гораздо приятнее, чем первое - я потратила на него всё утро и весь вечер воскресенья, но (надеюсь) что-то адекватное получилось, и почти всё даже более или менее понятно) Изначально очень запутывало употребление букв х, y и w всех по отношению к функции от w - может, лучше бы их было по-другому называть.

По датасету - было бы круто использовать какие-то более биологические данные, не про квартиры, а условно про транскриптомы или длины параподий - потому что покемоны это здорово, но пока непонятно, как именно это можно использовать)